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Abstract 

We present the numerical solution of the non-linear evolution equation for DIS on 
nuclei for x = 10"'^ 10""^. We demonstrate that the solution to the non-linear evolution 
equation is quite different from the Glauber - Mueller formula which was used as the initial 
condition for the equation. 

We illustrate the energy profit for performing DIS experiments on nuclei. However, it 
turns out that the gain is quite modest: xau — Sa^proton for the same parton density. 

We find that the saturation scale oc Aa . For gold the saturation scale Qs,Au ~ 
1.5 GeV at X = 10^'^. Such a large value leads to considerable contribution of the high 
density QCD phase to RHIC data and reveals itself in essential damping for both xGa and 

F2A. 
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1 Introduction 



Deep inelastic scattering on nuclei gives us a new possibility to reach a high density QCD phase 
without requiring extremely low values of x (extremely high energies). This can be seen directly 
from simple estimates of the packing factors for the partons in DIS on nuclei. Indeed, the packing 
factor is equal to 
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« A-^KNix,Q') (1.1) 

For a gold nucleus the packing factor is the same as for a nucleon at xa ~ 200 ■ x^, assuming 
xGNix-iQ"^) OC {1/x)^ with A ~ 0.3 which follows from the HERA experimental data [||]. The 
gluon density xG^^^^^ is the solution of the linear DGLAP evolution equation 0. The equation 
(|1.1| ) implies another advantage of nuclear targets. Namely, the value of the saturation scale is 
larger for nuclear targets than for nucleon ones. Indeed, the simplest estimate for the saturation 
scale can be deduced from the equation ka{x,Q1 ^i.^)) = 1- The solution of this equation is 

Ql^x) OC AW^ , (1.2) 

where the anomalous dimension 7 = dhi{xG isi) / dlnQ"^ is a smooth function of x and Q^. The 
fact that the saturation scale is larger for nuclei from the theoretical point of view makes our 
calculations more reliable. 

Determination of the saturation scale Qs,a{x) is one of the most challenging problems in QCD. 
The DGLAP equation [Q] describes the gluon radiation which increases the number of partons. 
However, when the parton density becomes large the annihilation processes enter the game and 
they suppress the gluon radiation. These effects tame the rapid increase of the parton densities at 
the saturation scale Qs,a{.x) 0, H, ||. At this scale the nonlinear effects become important and the 
parton evolution cannot be described by a linear equation any more. This argument stimulated 
a development of new methods @, |, |, |, 0, |, |, |10|, [11], 0, 0, |15|, g , which finally lead to the 



very same nonlinear evolution equation for the imaginary part of the DIS elastic amplitude^: 
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A^a(xoi,F;&) = A^A(xoi,>^o;&)exp 
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^ Eq. (^]^) was proposed in momentum representation by Gribov, Levin and Ryskin g and it was proved in 
double log approximation of perturbative QCD by Mueller and Qiu Q|, in Wilson Loop Operator Expansion at 
high energies by Balitsky in colour dipole approach jl^ to high energy scattering in QCD by Kovchegov (Oj, 
in effective Lagrangian approach for high parton density QCD by lancu, Leonidov and McLerran (see Ref. ]14| ] 
and Refs. for previous efforts). Braun derived this equation by summing "fan" diagrams, as GLR did, but 
using the triple BFKL ladder vertex of Refs. [Q. Therefore, this equation is a reliable tool for an extrapolation 
of the parton distributions to the region of low x. 
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The equation is written for N ^iju ^\ b) = Im adipoie('^-L! ^5 where Odipoie the amphtude of the 
elastic scattering for the dipole of the size r^. The total dipole cross section is given by 



(^d\po\e{r±,x) = 2 y d b NA{r±,x;b) . 
The total deep inelastic cross section is related to the dipole cross section 
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where the QED wave functions \E''^* of the virtual photon are well known |jT2|, |T8|, |T9l . The equation 
( |1.5| ) is quite transparent and it describes two stages of DIS The first one is the decay of 
a virtual photon into colourless dipole {qq -pair) which is described by the wave function in 
Eq. (11.51). The second stage is the interaction of the dipole with the target (cr^poie Eq. (|1.5D). 
This equation is the simplest manifestation of the fact that the correct degrees of freedom at high 
energies in QCD are colour dipoles |T^ . 

In the equation (O), the rapidity Y = — Inx and Yq = —Iuxq. The ultraviolet cutoff p is 
needed to regularize the integral, but it does not appear in physical quantities. In the large Nc 
(number of colours) limit Cp = Eq. ( p..3| ) has a very simple meaning: the dipole of size 

xio decays into two dipoles of sizes X12 and X02 with the decay probability given by the wave 

function l^'P = -t-^t-. These two dipoles then interact with the target. The non- linear term 
''02 ■"■12 

takes into account the Glauber corrections for such an interaction. 

The equation ( |1.3D does not depend on the target explicitly. This independence is a direct 
indication that the equation is correct for any target in the regime of high parton density. The 
only dependence on a target comes from initial conditions at some initial value Xq. For a nucleus 
target it was proven in Refs. |]^, 13] that the initial conditions should be taken in the Glauber 
form: 



A^a(xoi,xo;&) = A^GA/(xoi,a;o;&) 



with 



A^^Af(xoi,x;6) = 1 - exp 
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The cross section of the dipole interaction with a nucleon is defined as 
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with that has been calculated in the Ref. [^. The equation (|1.7|) represents the Glauber 
- Mueller (GM) formula which accounts for the multiple dipole-target interaction in the eikonal 
approximation ||18|, p2|, p3|. 
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It turns out, however, that at sufficiently high initial Xq {xq ~ 10 ^) the equation ( |1.7|) leads to 
the same initial profile as the simplified version 

iV^A^(xoi,x;6) = 1 - exp 



'^5^ ^01 ^qDGLAP 
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The nucleon gluon distribution inside a nucleus xG^*^^^^ is parameterized according to LO 
solution of the DGLAP equation 0. The function SA{b) is a dipole profile function inside the 
nucleus with the atomic number A. The value of xq is chosen within the interval 

exp(- — ) < xo < -\- , (1.9) 
as zmuA 

where Ra is the radius of the target. In this region the value of Xq is small enough to use the 
low X approximation but the production of the gluons (colour dipoles) is still suppressed since 
a5ln(l/2;) < 1. Therefore in this region we have instantaneous exchange of the classical gluon 
fields. Due to this fact an incoming colour dipole interacts separately with each nucleon in a 
nucleus 0. 

Solutions to the equation ( |T73D were studied in asymptotic limits in the Refs. p|, [l^, A 



first attempt to estimate numerically the value of the shadowing corrections in DIS on nuclei in 
the framework of the non-linear evolution equation was made in the Ref. 0. In that paper the 
AGL evolution equation was solved in the semi-classical approximation. It turns out that the 
AGL equation is the same as Eq. (|1.3|) in the double log approximation but it is written not for 
the amplitude N but for the opacity Q {N = 1 — exp(— f2/2)). Braun solved numerically an 



equation similar to Eq. (|1.3| ) in momentum representation with some specific initial conditions 
and we will discuss his results below. In our recent paper we solved Eq. ( |1.3| ) numerically for 
the proton targets. 



Our approach to the solution of Eq. ( |1.3|) is based on space representation. We have several 



arguments for working in this representation. First of all, it allows us to formulate the initial 
conditions. Moreover, the initial conditions are not known at large distances at all. Consequently 
we would expect problems calculating the amplitude in the momentum representation since it 
implies knowledge of the amplitude in the whole kinematic region of rj_. Our last argument is 
that in the space representation the unitarity constraints have the simplest form: 

2N{r^,x;b) = N\r^,x;h) + G,„(rx,x;6), (1.10) 

where Gj„ stands for contributions of all inelastic processes. We assume in Eq. ( |1.10| ) that at 
high energies the elastic amplitude is dominantly imaginary a'^'(r_L, x; 6) — > N{r±,x;b). The 
equation (|1.10|) implies N{x,r±;b) < 1 with a natural limit N{x,r±;b) — 1 at x — * 0. This 
limit provides us with a very useful check in our numerical approach. 

In the present paper we report on our solution of the equation (|1.3|) for the nuclear targets. The 



method of iterations proposed in is applied. Contrary to the Ref. |T^ we study solutions for 



phenomenologically reasonable kinematic regions only, and for real nuclei. The main goal of the 
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research is determination of the saturation scale Qs,a{x) and its dependence on the target atomic 
number A. We also present estimates of the gluon density xGa and the structure function F2A 
in the region of the high parton density QCD. 

The paper is organized as follows. In the next section (2) we give a review of the nuclear target 
properties. In section 3 a solution of the nonlinear equation is presented. Section 4 is devoted to 
the calculations of xGa and F2a- In section 5 we determine the saturation scales and study their 
A-dependence. Section 6 presents discussion of the self-consistency of the numerical procedure 
and results. We conclude in the last section (7). 



2 Nuclear input data 



In this section we discuss the properties of the nuclei our computations are made for. The nuclei 
are described by their radial nucleon distribution p{r). This distribution is normalized to the 
total nucleon number A: 

J p{r)(fr = A. (2.11) 

The nucleon distribution for all nuclei discussed in this paper are parameterized by one of the 
following three parameterizations: 

1. Two-parameter Fermi model (2pF) 

P(^) = Po/ (1 + exp ((r - c)/z)) . 

2. Three-parameter Fermi model (3pF) 

p(r) = Po (1 + w r'^ / c^) / {1 + exp ((r — c)/z)) . 

3. Three-parameter Gaussian model (3pG) 

p{r) = po(l + wrVc^)/(l + exp {{r^ - c^) / z'^)) . 



The constants c, z, and w are free parameters varying from nucleus to nucleus. The table (|I|) lists 
the relevant input data for all nuclei used in this work p6|. The parameter po is always found 



due to the normalization ( 2.11 ). 

The radius r can be decomposed to the transverse b and longitudinal ri components 

2 1,2 I 2 
r =0 -I- r, . 

Then the transverse profile function 5*^(6) is defined as the integral of the radial density p over 
the longitudinal component r/: 

SA{b) = I p{r)dri. (2.12) 



Fig. I] displays the function 5*^(6) for various nuclei. 
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Ca 


40 


3pF 


3.486 


3 6685 


5839 


-0.1017 


Zn 


70 


2pF 


4.044 


4.409 


0.583 




Mo 


100 


3pG 


4.461 


4.559 


2.6723 


0.339 


Nd 


150 


2pF 


5.048 


5.7185 


0.651 




Au 


197 


2pF 


5.33 


6.38 


0.535 





Table 1: Parameters for the nucleon distributions in the nuclei taken from the Ref. [p^ ]. 
SA(b) 



Figure 1: The nucleon transverse profile function inside the 
nucleus Syi(6) is plotted function of b. 
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We wish to draw reader's attention to the fact that S{h = 0) for proton is almost the same as for 
A^e (see Fig. 0). This fact has an obvious explanation: partons inside the proton are distributed in 
smaller area with R ^ 3 GeV~^ = 0.6 fm which is smaller than the area given by extrapolation 
of the Wood-Saxon parameterization to A = 1. Indeed, Ra=i — 1.2 fm > R = 0.6 /m. This 
simple fact makes our estimates based on Eq. ( p..l|) too optimistic. We will show below that the 
packing factor for gold is the same as for nucleon at xau ~ ( 5 -i- 10 ) ■ xat . 



3 Solution of the nonlinear equation 



In the Ref. [^] we suggested to solve the equation (|1.3|) by the method of iterations. All details 
about our method can be found in that paper, where the equation (|1.3|) was solved numerically 
for the proton target. 



Following the very same procedure as in we obtain numerical solutions for all nuclei listed 
above. For xG^'^^"^^{x,Q'^) we use the GRV'94 parameterization and the leading order solution 
of the DGLAP evolution equation p^. The initial conditions (|1.6| ) are set at Xq = 10~^. The 
constant value for the strong coupling constant as = 0.25 is always used. The solutions are 
computed within the kinematic region down to x = 10"'' and distances up to a few fermi. 

The function A^^ is formally a function of three variables: the energy variable x, the transverse 
distance r±, and the impact parameter b. The 6-dependence is parametric only because the 
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evolution kernel does not depend on b. In order to simplify the problem we assume that the 
function preserves the very same 6-dependence as introduced in the initial conditions:^ 



-KA{x,r^)SA{b)/SAiO) 



)• 



(3.13) 



NA{r±,x;b) = (1 - 

The function ka is related to the "6 = 0" solution iV^(r_L,x): 

KA(x,r^) = - ln(l - NA{r^,x)). (3.14) 

NA{r±,x) represents a solution of the very same equation (|1.3| ) but with no dependence on the 
third variable. The initial conditions for the function NA{f±,x) are taken at 6 = 0. For the case 
of the proton target |^ the anzatz in the form (|3.13|) is shown to be a quite good approximation 
of the exact 6-dependence of the solution to ( |1.3| ). We will investigate the accuracy of this anzatz 
for gold in the end of this section. 

For each nucleus the function NA{r±,x) is obtained after about ten iterations. The Fig. ^ shows 
the solutions for various nuclei. Note that though the solutions are of similar form they differ 
from each other indicating some memory of the initial conditions. 
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Figure 2: The function Na is plotted versus distance. The four curves show the result for the nuclei 
Ne, Ca, Mo, and Au. 



It is worth to compare the obtained solutions with the Glauber formula ( |1.8| ). For example the 
comparison is made for two nuclei: the lightest one A^e and the most heavy Au. The two models 
are plotted together in the figure |^. 

^Note that for the Gaussian parameterization of the profile function 5*^, A = R\ Sa{Q). 
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Figure 3: The functions Na for Au and Ne (dashed lines) are plotted together with the corresponding 
Glauber formulas (dotted lines). 



Having obtained the solutions for nuclei we can answer the question about the energetic gain of 
performing DIS on nuclei. To this goal the results for nuclei should be compared with ones for 
proton targets ||2l| . Unfortunately, the gain is quite modest and it consists of about half order of 



magnitude for very heavy nuclei. Namely, for a gold target Au the scattering amplitude Nau (and 
consequently the packing factor kau) almost coincides with the proton one when xau — Sxproton, 



which is quite different from the previous optimistic estimates (see Eq. (IJ.) and discussion there). 
The simple reason for such a mild gain was mentioned in the previous section and it is related to 
the fact that Sat (6 = 0) is almost the same as for A^e as can be seen in the Fig. |1|. 

It is possible to determine the gluon density xGa from the solutions obtained. To this goal the 
accuracy of the anzatz ( ^.131 ) has to be estimated. The 6-dependence so far ignored in the function 
Na should be investigated. In order to restore the 6-dependence of the function Na we solve the 



equation (|1.3|) with the only assumption that the impact parameter b is much larger than the 
dipole sizes: 

xoi < b; Xo2 < b. (3.15) 

In this case the 6-dependence of the rhs of ( p..3|) is significantly simplified. We perform the 
computations for only one nucleus Au, which is the most different from the proton among the 
nuclei that we have considered. The Fig. ^ shows examples of the comparison between the correct 
6-dependence of the solution and the anzatz (Eq. ( ^.13| )). It can be seen from the figure Fig. ^ 
that our anzatz ( p.l3|) is a very good approximation for relatively small b and up to 6 ~ Ra- 



For larger impact parameters our anzatz underestimates the true solution. This fact can be 
simply understood theoretically. At very large impact parameters the exponential in ( |3.13|) can 
be expanded and we get Na = HAS{b)/ S{0). Then substitution of the anzatz to the equation ( pT3D 
shows that it underestimates the solution indeed. For such values of the impact parameter it is 
more natural to search for the true solution of ( |1.3|) in the form NA{r±, x; b) = n^(r_L, x) S{b) / S{0) 
and to write down an evolution equation directly for the function ha- In this case, we would obtain 



the function ua somewhat larger than our ka- Such an approach was adopted in the Ref. [|28| but 
for the whole range in b which seems to be wrong for small impact parameters and the Glauber 
initial conditions (O 
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Figure 4: The function Nau is plotted versus the impact parameter b. The solid line is the exact b 



dependence, while the dashed line is the anzatz (3.15). 
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4 Gluon density and the F2 structure function 



4.1 Gluon density xGa 



We can proceed now with the computation of the gluon density xGa which we define using the 
Mueller formula [0: 



oo (1^.2 

4/Q2 



(4.16) 



There is some uncertainty in ( |4.16| ) due to 6-integration. Supposing that the main contributions 
to the 6-integration come from the region b ~ Ra, the uncertainty can be estimated. In this region 
the anzatz (|3.13|) underestimates the solution by maximum about 25%, and hence the obtained 
gluon density is also expected to be underestimated by the same amount. The Fig. |^ presents 
results of the calculations for the gluon density xGa for three values of while the Fig. ^ shows 
its dependence for various x. 
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Figure 5: The gluon density xGa is plotted as a function o/lgx for the nuclei Ne, Ca, Zn, Mo, Nd, 
and Au. 
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Figure 6: The gluon density xGa is plotted as a function of Q^; 
X = 10-^; d- x = 10~3. 
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Table 2: The power n for various values 
of X. 



The gluon density xGa defined as in ( [4.16|) grows with decreasing x. In a sense, this observation 
contradicts the "super-saturation" of Ref. |T6[ where the density is predicted to vanish in the 
large rapidity limit. However, the gluon density definition accepted in Ref. [jl6l differs from ours, 
and hence the functions are not compatible. 

It is important to investigate the A dependence of the gluon density xGa- The leading twist 
perturbative QCD prediction based on DGLAP equation is xGa = AxGn, where xGn is a 
nucleon gluon density. 

Let parameterize the A dependence as xGa = A^-xG^. The results obtained for the power n are 
presented in the table (0). At high x we indeed obtain the law xGa — AxG^ which is shown 
in the Fig. |[ At moderate x ~ 10"^ a transition occurs and n turns out to be quite different 
from 1. At low x the gluon density xGa is rather proportional to A"^^^ (Fig. For large the 
transition occurs at somewhat smaller x, but the dependence on is very weak (Table @). This 
described A dependence is actually anticipated from general theoretical arguments (see Ref. ^ 
and references therein). 



4.2 Anomalous dimension 

We define the average anomalous dimension as follows 



7{A,x,Q') 



2, d(\n xG a{x , Q"^)) 



rflnQ2 



(4.17) 
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The results of the computations are presented on the Fig. ^a. To our surprise the anomalous 
dimension almost does not depend on the atomic number A displaying independence on the target. 
However, the anomalous dimension for the function satisfying Eq. (|1.3| ) shows a considerable 
A dependence as can be seen in the Fig. |^ (b-c). We define 



In 



d\nNAi2/Q,x) 
dlnQ2 



+ 1 



Such a different behavior of 7 and 77V, in our opinion, is mostly related to the fact that shorter 
distances enter the calculation of 7 compared to 7Ar as can be observed directly from Eq. ( |4.17] ). 
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Figure 7: The anomalous dimensions are plotted as a function of Igx. (a) - j for the various values of 
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4.3 F2A 

The structure function F2A is related to the dipole cross section 



V POO (ir f 

F2Aix,Q') = - — d%NAir,x;b). 



(4.1^ 



The Fig. ^ presents the function F2A for various nuclei. 

We now in the position when we can try to compare our results with ones obtained in the Ref. 
|T^. That paper presents results for the structure function F2 computed for the lead nuclei (A 
= 207). Unfortunately, the kinematic regions investigated in our work and in the Ref. [|16| have 
a very small overlap. Nevertheless at not extremely low x and not too high < 1000 GeV^ a 
comparison can be made (we use our gold [A = 197) calculations). We fail to reproduce results 
of the Ref. 



16| . In all the region of the comparison we predict a few times smaller values for the 



function F2a- 
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Figure 8: The structure function F^aIA is plotted as a function of Ig x for the nuclei Au, Ala, and Ne. 
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There are several reasons for the obtained mismatch. First of all, there are few minor ones 
related to the difference between the nuclei compared and different as used in the computations. 
In addition, our treatments of the 6-dependence are quite different. However, in our opinion, the 
main reason is in the initial conditions of the evolution, which do not coincide with the ones of 



the Ref. |]16|. At moderate a; (x = 10 ^10 ) our solution depends on initial condition as well 
as at lower x where such a dependence is concentrated at very short distances. 

Concluding this comparison we would like to make the following comment. We actually believe 
that our results on F2A are more reasonable. At not very low x ~ 10~^ ^ 10~^ we expect 
the function F2A to be proportional to A. Namely, F2A — AF2. Our results on F2A display 
this behavior and are in a good agreement with the experimental data on the proton structure 



function F2. Contrary to us, the results of [jT6[ seem to be few times off from this agreement. 



5 Saturation 

It is natural to define the saturation scale through the equation 

ka{x,2/Qs) = 1/2. (5.19) 
The definition ( 5.19| ) agrees with the one adopted for the Glauber formula (|1.8|) . 



Fig. ^ displays the saturation scale Qs,a obtained in ( p.l9[ ) for various nuclei. 

The main question which we want to study is the dependence of the saturation scale Qs,a{x) on 
the atomic number A. Let us assume a power law behavior for the saturation scale Qs,Aix) as a 
function of A: 

QsA^) = Cix)AP^^\ (5.20) 

where C{x) and p{x) are x dependent functions. The power p{x) is of our main interest. In order 
to check the anzatz ( ^.201) and find the power we study the A dependence of the saturation scale 
in double logarithmic scale and find it to fit well a straight line that justifies use of the anzatz 
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Figure 9: The saturation scale Qs,A "is plotted versus Igx (a) and Ra (b). 



Nuclei \ X 


10-^ 


10-^ 


10-^ 


10-^ 


10-3 


Light 


0.18 


0.22 


0.26 


0.31 


0.49 


Heavy 


0.15 


0.19 


0.22 


0.24 


0.31 


All 


0.17 


0.20 


0.24 


0.27 


0.39 



Table 3: The power p{x) for various val- 
ues of X. 



(|5.20| ). The power p{x) can be found by the least square fit. The table presents the fit for the 
three cases: light nuclei (A^e, Ca, Zn); heavy nuclei {Zn, Mo, Nd, and Au); all nuclei together. 

The table (0) is quite transparent. The power p{x) decreases with decreasing x. At the beginning 
of the evolution the A-dependence of the saturation scale for the light nuclei is close to the law 
Qs,A ~ v4^/^, at X ~ 10~^ the saturation scale Qs,a ~ A^^^, while at higher energies it tends to 
Ql A ^ A^^^- For the heavy nuclei the situation is similar but the decrease of the power p(x) 
is significantly slower. The above observations are in complete contradiction with conclusions 



derived in the Ref . ||2^ , where the saturation scales were deduced from the equation ( |1.3| ) in the 



double logarithmic approximation. The main source of this large discrepancy is the fact that the 
anomalous dimension in the solution of the DGLAP equation turns out to be larger than 1/2 
which is the maximum value for the BFKL evolution in the leading order (see Eq. ( |1.2| )). This is 
a disturbing result since it could mean that the leading order non-linear evolution equation (see 
Eq. (|1.3| )) is not enough to make a reliable predictions in the kinematic region of high density 
QCD. 

We would like to stress that the numbers presented in the table (H) are quite approximate. These 
numbers display a certain sensitivity to the definition of the saturation scale (we use Eq. ( p.l9| ) 
but other definitions can be investigated as well). Moreover, the powers obtained are results of 
regressions over too small number of points, which actually implies large errors. However, we are 
convinced that the decreasing property of the power p{x) with high energies is quite general. 
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Table @ can be explained qualitatively in terms of the anomalous dimension 7. The following 
arguments are presented for the GM formula (Eq. ( |1.8[ )) or /and for the packing factor (Eq. ( |1 .11 )), 
but the physical picture is valid for the solution of the equation ( |1.3|) as well. In the GM formula 
the power p(x) = 0(13:^ • Though 7 is usually assumed to be small, it is often not the case. 
Moreover, 7 is actually x and Q dependent. At large x and correspondingly small saturation 
scale Qs,A the GRV gluon density parameterization implies very large 7 tending to unity (see e.g. 
Ref. 0). This is the origin of the large powers in the table (0). When we go from light nuclei 
to more heavy (at fixed x) the saturation scale obtained increases, and consequently 7 decreases. 
That is why smaller powers are obtained for the heavy nuclei compared to the light ones. As x 
decreases, on one hand 7 increases at fixed Q. On the other hand, at smaller x a larger saturation 
scale Qs,A is obtained. Finally, with decreasing x the power p{x) decreases as well tending to the 
value p{x) = |. 



6 Discussion of the results 

Let us explain qualitatively the A dependence of the results. We also wish to perform checks for 
the self-consistency of the results obtained. 

The main starting observation is that the anomalous dimension 7 ( |4.17|) is almost A independent 
while its dependence become weaker as x decreases. 

Define the power a: 



.(A:..Q') = - '""';;t"^y^» . (6.21) 



We find that a ~ 1 — 7. Consequently we can try to write ka in the power-like form: 

fiA{x,r^) ~ A^(^)(r2)". (6.22) 
In the equation (|6.22|) we introduced the A-dependence in order to define the saturation scale: 



QsA^) ~ ^^^^'"^ = ^"^"^ • (6-23) 

The power (3 is obtained to be almost independent. It ranges from about 0.33 at x = 10~^ and 
to about 0.17 at a; = 10~^. Meanwhile at the saturation scale a is almost constant a ^ 0.45-4-0.5. 
This corresponds to 7 ^ 0.5, which is the anomalous dimension of the linear BFKL term in the 
nonlinear equation ( |1.3| ). 

Assuming be relatively small we would obtain for the gluon density 



xGa{x,Q^) ^ R\ — —KA{x\r). (6.24) 

Jx X' JA/Q'^ 
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Substituting (|6.22| ) to the equation ( |6.24| ) we obtain 



xGa{x, Q') ^ R\ {Q^f-- . (6.25) 

J X X 



Assuming that the integral in ( |6.25| ) is dominated by x' ~ x, we see that within the above 



approximation 7 simeq 1 — a indeed. The independence of the anomalous dimension on A is a 
consequence of the fact that P happened to be almost independent. 

What about A-dependence of the gluon density? We obtain for the power n the relation 

n{x) = 2/3 + P{x) = 2/3 + 2p{x)a{x) . (6.26) 

The equation ( |6.26D is in a sense an additional consistency check. The numbers presented in 
the tables (|[) and (0) are in a good agreement with each other (though the equation ( |6.26|) is 
not fulfilled exactly due to many approximations done on the way to ( |6.26| )). The fact that the 
powers n in the table (|^) do not much sensitive to is now understood as a consequence of the 
corresponding independence of /3. 



7 Conclusions 

In this paper we reported on the exact numerical solutions of the nonlinear evolution equation 
(|1.3| ) for nuclear targets. The solutions are obtained by the method of iterations proposed in 



the Ref. [^. Various nuclei starting from the very hght Ne2o and up to heavy Auigy were 
investigated. 

We demonstrated that the solution to the non-linear equation is quite different from the GM 
model that has been used for estimates of the saturation effects. However, this model can be used 
as a first iteration of the non-linear equation which leads to faster convergence of the numerical 
procedure. 

From the experimental point of view the obtained results support the energetic profit for per- 
forming DIS experiments on nuclei. However, the gain is quite modest, and it can be estimated 
about five times for gold targets compared with proton. 

The gluon density xGa and the structure function F2A were estimated. At small x ^ 10^^ ^ 10^^ 
the damping due to non-linear effects leads to suppression of a factor 2^4 for heavy nuclei in 
comparison with the DGLAP prediction xGa = AxG^'^^^^ . At moderate x ~ 10~^ -r- 10~^ the 
structure function F2A — AF2, where the obtained values for F2 agree well with the experimental 
data on proton. 

The dependence of the gluon density xGa was investigated as a function of the atomic number 
A. At high X, the density is shown to be proportional to A. A transition occurs at moderate 
X ~ 10~^. At small x, xGa is rather proportional to A^^^. To our knowledge this is a first 
numerical confirmation of this dependence from the master equation Eq. (|1.3|) expected on general 
grounds. 
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The found solutions of the nonhnear equation were used to estimate the saturation scale Qs,a{x)- 
In agreement with all theoretical predictions the saturation scale grows with decreasing x. 

The dependence of the saturation scales on the atomic number A was a focus of the research. A fit 
of the saturation scale to the power law Qs,a{x) ~ A^^^^^ was investigated. The numerical values 
obtained for the power p{x) are actually sensitive to the saturation scale definition. Nevertheless, 
we predict a decreasing behavior of the power with decreasing x for both heavy and light nuclei. 

One of the interesting properties of the solutions of the equation ( p..3|) is that they display the 
scaling phenomena. Namely, the solution A^^ is not a function of two variables x and r± but 
rather a function of a single variable r±Qs^A{x)- A complete report on this subject is now in a 
final stage of preparation and will appear shortly ||30[| . 

Among several observations which come from the numerical solutions we would like to point 
out that the weak dependence of xGa anomalous dimension on the atomic numbers as well as 
the fact that ^ (x A^^^ at a; — ^ look quite impressive for us. They certainly need some 
theoretical explanation and we hope they will stimulate such explanations. We firmly believe 
that our calculations of Qs,a for various nuclei will provide a theoretical basis for the discussion 
of the RHIC data in high parton density QCD approach . 



When the present paper was finished for publication we read a new paper devoted to the very 
same subject. We have to stress that the initial conditions of our analyses are quite different. 
Though within the evolution in x the influence of initial conditions become weaker they are 
important for quantitative analysis in the experimentally reasonable region of not too small x. 
We were pleased to discover that despite the difference in the initial conditions the saturation 
scale Qs obtained in PT| is proportional to A^^^, which is a dependence quite similar to ours. 



The scaling phenomena confirmed by the numerical calculations of the Ref. has been also 
observed by one of us and the report will appear soon [Q. 
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